ICFEP/ABAQUS-UQ

Consisting of:

Enriched experimental of design + Bayesian inference

clearvars; clc; close all;

Part 1: Setup

Directory setup

%% Add UQlab core and subroutine directory
addpath("E:\OneDrive - Imperial College London\PhD_Work\UQlab\UQLab_Rel2.1.0\UQLab_Rel2.1.0\core");%UQLAB core
addpath("E:\OneDrive - Imperial College London\PhD_Work\SBI\Main\subroutines");%subroutine
 
% tell whether in the running folder
if ~contains(pwd(),caseInsensitivePattern('run'))
error('Working path is wrong!!!');
end
 
 
% define a working directory structrure for input/output
% decide whether reserve current FE results: true <-- delete; false <-- reserve
directory = Def_Dirs('true');
 
%% Initialize UQlab (not printing in the command window)
seed = 4;
rng(seed,'twister');
evalc('uqlab;');
 

Probabilistic input for EoD/parameters/discrepancy setup

%% Read input EoD parameters and observation error discrepancy hyperparameters range from csv
[Prior_forward,Prior_hyperdiscrep,Prior_EoD,AnParam] = readPrior(directory);
Parameter priors for forward predictive purpose
Prior_forward = 3×1 struct
FieldsInput_nameparameter1parameter2type
1K_{0}1.351.8uniform
2OCR1550uniform
3G_{0}55000165000uniform
Parameter priors for experimental design
Prior_EoD = 3×1 struct
FieldsInput_nameparameter1parameter2type
1K_{0}1.351.8uniform
2OCR1550uniform
3G_{0}55000165000uniform
Hyper parameter (Gaussian-type observation error) for outputfield 1
ans = 5×4 table
 Input_nameparameter1parameter2type
1"w0"010"uniform"
2"w1"-1010"uniform"
3"w2"-1010"uniform"
4"w3"-1010"uniform"
5"h"020"uniform"
Hyper parameter (Gaussian-type observation error) for outputfield 2
ans = 1×4 table
 Input_nameparameter1parameter2type
1"w0"010"uniform"

Observation data setup

% if this is a systhetic problem, add Gaussian noise percentage
Noise.scalar = 0;% Value as a percentage (e.g. 0.01-1.0)
 
%% read observation (ground truth) data
[observationcoord,observationsarray,gtarray,Noise,observation_stageTime,AnParam] = readobservation(directory,Noise,AnParam);

Analysis parameters setup

AnParam.FE_package = 'ICFEP';% choose FE package
AnParam.EoDseed = seed; %restore the seed
AnParam.N_RUN_EoD = 30; %number of FE runs required for EoD
AnParam.N_FE_stage = 70; %number of FE loading stages
AnParam.N_enrich = 71; %hope to enrich N_enrich-th stage for surrogate construction
% (1≤ AnParam.N_enrich ≤
% AnParam.N_FE_stage, if not in this range, no enrich)
AnParam.Switch_Predictive = "on"; % switch for prior/posterior predictive for real outputs
AnParam.TestDataRun = 5; %Test FE runs number
AnParam.TrainDataPerc = 1; %rangePerc = (e.g., 0.8-1), TrainRun = rangePerc * (Allrun - TestDataRun);
AnParam.Switch_surrogatePlot = "on"; % switch to show the test error/reconstruction error (only show in mlx)
AnParam.Switch_corner_scatter = "on"; % switch to show the posterior corner figures
AnParam.Switch_Sobol = "on"; % switch to show the posterior corner figures
AnParam.Switch_threeSigma = "on"; % switch to show the posterior paramters 3 sigma with stages
AnParam.Switch_RMSE = "1D"; % switch to show the posterior paramters mean and std with stages
AnParam.Trim = 1:size(observationcoord{1},2); %Trim range for FE output after interpolation
AnParam.DR = "on"; % choices for surrogates
AnParam.Surrogate = "PCE"; % choices for surrogates
AnParam.Switch_priorBound = "on"; % switch for prior forward boundary
AnParam.inferredCopula = "off"; % this button always encounter error with AnParam.Switch_priorBound
AnParam.figureDisplay = "Show-in-mlx"; % switch for display/save figs
AnParam.LL_Normalize = "off"; % normalize LL by MAP values
AnParam.LL_MAP_test = "off"; % test LL evaluation values at MAP and save in a csv file
 
 
 
% Used in readoutput.m Very difficult to get rid of this as .m is hidden by UQlab.
if strcmp(AnParam.FE_package,'ICFEP')
%%% Define global variables
global FE_ID
FE_ID = 0;
end

UQlink setup

%% Wrapper ICFEP/ABAQUS with UQLink:
modelname = 'NonLinear_Truss';
filename = 'NonLinear_Truss';
switch AnParam.FE_package
case 'ABAQUS'
myUQLinkModel = ABAQUS_wrapper(modelname, filename);
case 'ICFEP'
myUQLinkModel = ICFEP_wrapper(modelname, filename);
end

Part 2: EoD + Bayesian inference

2.2 Recursive inference

for jj = 1:AnParam.N_Obs_stage
 
% enrich EoD manually and extract enriched FE restults
if jj==1 %First stage EoD run
[~] = EnrichedFERun(myUQLinkModel,Prior_EoD, AnParam,filename,directory,jj);
else
[~] = EnrichedFERun(myUQLinkModel,Prior_EoD, AnParam,filename,directory,jj, Posterior);
end
 
% Merge different EoD results to EDarray/FEcell_output if possible
% Update the number of FE runs if possible
[coord, EDarray, FEcell_output,AnParam,FE_stageTime] = merge_EoD(AnParam);
 
 
%%Note: (1) FE_stageTime and observation_stageTime may be not consistant,
% (2) FE coordination and observation coordination may be not consistent
%%interpolation to make FE result to fit with observation
%%(FE_stageTime/FE coordination must contain observation_stageTime/observationcoord)
FEcell_output_interp = interp_coord_stageTime(observationcoord,coord,FEcell_output, ...
AnParam,observation_stageTime,FE_stageTime);
%% Trim coordination and FE results if possible (trim outputput 1)
[coord,FEcell_output_interp,observationsarray] = Trim(jj,AnParam,FEcell_output_interp, ...
coord,1,observationsarray);
 
%% Construct PCA-PCE for a given stage {jj-th stage}
[myPCE, PCA] = PCEPCA(AnParam, EDarray, FEcell_output_interp, jj);
 
%% Make PriorOpts/myPriorDist structure for parameter/observation discrepancy priors
%% Initial stage prior should be defined manually
 
if jj ==1
%intial srage prior should be defined manually
PriorOpts = []; Posterior{jj} = []; %initial empty cells
[PriorOpts, myPriorDist] = Dynamic_PriorOpts(Prior_forward,jj,AnParam, ...
observationsarray,Prior_hyperdiscrep,Posterior);
else
%rest stages should be inferred automatically
PriorOpts = []; %clear PriorOpts from last stage
[PriorOpts, myPriorDist] = Dynamic_PriorOpts(Prior_forward,jj,AnParam, ...
observationsarray,Prior_hyperdiscrep,Posterior);
end
 
%%create empty folderes to restore results at initial stage
[~] = emptyFolder(jj)
 
%%%prior predictive
[~] = predictive_prior(AnParam,EDarray,FEcell_output_interp,coord,observationsarray, ...
gtarray,observation_stageTime,directory,myPriorDist,jj);
 
%% === Observations ====
%%UQLAB doesnot support NAN in the myData.y
%%Repalce NAN with numerical value eplision=1e-20
epsilon = 1e-20;
for kk = 1:AnParam.N_outputfields
if ~isnan(observationsarray{kk}(:,jj))
obsconcat{kk} = observationsarray{kk}(:,jj)';
else
obsconcat{kk} = epsilon*ones(size(observationsarray{kk}(:,jj)))';
end
end
myData.y = obsconcat;
myData.Name = 'Multiple observations data sources';
 
 
% %% === Likelihood switch ====
myLogLikeli = likeliSwitch(jj,observationsarray,PCA,myPCE,AnParam,Prior_hyperdiscrep,coord,directory);
 
 
%Add boundary to forward prior (since UQlab not supporting trunction during coupula/marginal input
% inference) we need put additional boundary afterwards
myPriorDist = boundary(PriorOpts,Prior_EoD,AnParam,myPriorDist);
 
 
%% Solver options: AIES MCMC
Solver.Type = 'MCMC';
%Solver.MCMC.Visualize.Parameters = 1:AnParam.N_parameters;
%Solver.MCMC.Visualize.Interval = 20;
Solver.MCMC.Sampler = 'AIES';
Solver.MCMC.Steps = 300;
Solver.MCMC.NChains = 100;
BayesOpts.Data = myData;
BayesOpts.Type = 'inversion';
BayesOpts.LogLikelihood = myLogLikeli;
BayesOpts.Solver = Solver;
BayesOpts.Prior = myPriorDist;
BayesAnalysis = uq_createAnalysis(BayesOpts);
 
%% response surface for RSM
[~] = RMSE_plotting(BayesAnalysis,AnParam,jj,myPCE, ...
PCA,EDarray,observationsarray,gtarray,6000);
 
 
%% sobol indices for plotting
[~] = Sobol_indices_plotting(AnParam, FEcell_output_interp,EDarray,jj,observation_stageTime);
 
%% Postprocess (70% burn in; extract samples from posterior )
Posterior = BI_postprocess(jj,BayesAnalysis,AnParam,Posterior,gtarray,myPriorDist,Prior_forward);
 
%% save/display MAP for each stage with stages going
[~] = MAP(AnParam,Posterior, directory);
 
%% plot paramter three sigma range with stages going
[~] = threeSigma_plotting(AnParam,Posterior, gtarray,jj,Prior_forward);
 
%% Predictive posterior
[~] = predictive_posterior(Posterior,AnParam,jj,EDarray,FEcell_output_interp, ...
coord,observationsarray, gtarray,observation_stageTime,directory);
 
%% save LL evaluation values at MAP in a csv file if required
[~] = LL_MAP_test(BayesAnalysis,AnParam, jj)
 
 
end
PCA reconstruction MAPE is 0.5002%
PCA-PCE test MAPE is 1.9555%
PCA reconstruction MAPE is 0%
PCA-PCE test MAPE is 0.12902%
Starting AIES... |##############################| 100.00% Finished AIES!
Warning: When plotting RSM-response surface for outputfield 1, no observation for this stage!!
Print_MAP = 1×3 table
 K_{0}OCRG_{0}
11.350030.99851.0594e+05
Analysis for Stage 1 is finished! ╔═══╗───────╔╗──────╔╗ ║╔══╝───────║║──────║║ ║╚══╦╦═╗╔╦══╣╚═╦══╦═╝║ ║╔══╬╣╔╗╬╣══╣╔╗║║═╣╔╗║ ║║──║║║║║╠══║║║║║═╣╚╝║ ╚╝──╚╩╝╚╩╩══╩╝╚╩══╩══╝ ======================================================================================================================================================
PCA reconstruction MAPE is 1.1811%
PCA-PCE test MAPE is 6.4247%
PCA reconstruction MAPE is 0%
PCA-PCE test MAPE is 0.55001%
Warning: Maximum likelihood estimation did not converge. Iteration limit exceeded.
Warning: Maximum likelihood estimation did not converge. Iteration limit exceeded.
Warning: Unable to compute a covariance matrix because the computed Hessian matrix is not positive definite.
Warning: Maximum likelihood estimation did not converge. Iteration limit exceeded.
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.419306e-18.
Warning: Maximum likelihood estimation did not converge. Iteration limit exceeded.
Warning: Unable to compute a covariance matrix because the computed Hessian matrix is not positive definite.
Warning: cannot fit gumbelmin distribution to the given data
Warning: Maximum likelihood estimation did not converge. Iteration limit exceeded.
Warning: Maximum likelihood estimation did not converge. Iteration limit exceeded.
Warning: Unable to compute a covariance matrix because the computed Hessian matrix is not positive definite.
Warning: cannot fit gumbel distribution to the given data
Warning: cannot fit gumbelmin distribution to the given data
Starting AIES... |##############################| 100.00% Finished AIES!
Warning: When plotting RSM-response surface for outputfield 1, no observation for this stage!!
Print_MAP = 2×3 table
 K_{0}OCRG_{0}
11.350030.99851.0594e+05
21.373543.65851.0414e+05
Analysis for Stage 2 is finished! ╔═══╗───────╔╗──────╔╗ ║╔══╝───────║║──────║║ ║╚══╦╦═╗╔╦══╣╚═╦══╦═╝║ ║╔══╬╣╔╗╬╣══╣╔╗║║═╣╔╗║ ║║──║║║║║╠══║║║║║═╣╚╝║ ╚╝──╚╩╝╚╩╩══╩╝╚╩══╩══╝ ======================================================================================================================================================
PCA reconstruction MAPE is 1.1945%
PCA-PCE test MAPE is 4.0286%
PCA reconstruction MAPE is 0%
PCA-PCE test MAPE is 0.65542%
Warning: Maximum likelihood estimation did not converge. Iteration limit exceeded.
Warning: Unable to compute a covariance matrix because the computed Hessian matrix is not positive definite.
Warning: Maximum likelihood estimation did not converge. Iteration limit exceeded.
Warning: Unable to compute a covariance matrix because the computed Hessian matrix is not positive definite.
Warning: cannot fit student distribution to the given data
Warning: Maximum likelihood estimation did not converge. Iteration limit exceeded.
Warning: Unable to compute a covariance matrix because the computed Hessian matrix is not positive definite.
Starting AIES... |##############################| 100.00% Finished AIES!
Warning: When plotting RSM-response surface for outputfield 1, no observation for this stage!!
Print_MAP = 3×3 table
 K_{0}OCRG_{0}
11.350030.99851.0594e+05
21.373543.65851.0414e+05
31.404120.43716.7718e+04
Analysis for Stage 3 is finished! ╔═══╗───────╔╗──────╔╗ ║╔══╝───────║║──────║║ ║╚══╦╦═╗╔╦══╣╚═╦══╦═╝║ ║╔══╬╣╔╗╬╣══╣╔╗║║═╣╔╗║ ║║──║║║║║╠══║║║║║═╣╚╝║ ╚╝──╚╩╝╚╩╩══╩╝╚╩══╩══╝ ======================================================================================================================================================
PCA reconstruction MAPE is 1.1014%
PCA-PCE test MAPE is 1.9846%
PCA reconstruction MAPE is 0%
PCA-PCE test MAPE is 0.70975%
Warning: Maximum likelihood estimation did not converge. Iteration limit exceeded.
Warning: cannot fit gumbel distribution to the given data
Warning: cannot fit gumbelmin distribution to the given data
Warning: Maximum likelihood estimation did not converge. Iteration limit exceeded.
Warning: Maximum likelihood estimation did not converge. Iteration limit exceeded.
Starting AIES... |##############################| 100.00% Finished AIES!
Warning: When plotting RSM-response surface for outputfield 1, no observation for this stage!!
Print_MAP = 4×3 table
 K_{0}OCRG_{0}
11.350030.99851.0594e+05
21.373543.65851.0414e+05
31.404120.43716.7718e+04
41.675622.18018.3515e+04
Analysis for Stage 4 is finished! ╔═══╗───────╔╗──────╔╗ ║╔══╝───────║║──────║║ ║╚══╦╦═╗╔╦══╣╚═╦══╦═╝║ ║╔══╬╣╔╗╬╣══╣╔╗║║═╣╔╗║ ║║──║║║║║╠══║║║║║═╣╚╝║ ╚╝──╚╩╝╚╩╩══╩╝╚╩══╩══╝ ======================================================================================================================================================
PCA reconstruction MAPE is 4.6992%
PCA-PCE test MAPE is 1.8214%
PCA reconstruction MAPE is 0%
PCA-PCE test MAPE is 0.66844%
Warning: cannot fit gumbel distribution to the given data
Warning: cannot fit gumbelmin distribution to the given data
Warning: cannot fit gumbel distribution to the given data
Warning: cannot fit gumbelmin distribution to the given data
Starting AIES... |##############################| 100.00% Finished AIES!
Print_MAP = 5×3 table
 K_{0}OCRG_{0}
11.350030.99851.0594e+05
21.373543.65851.0414e+05
31.404120.43716.7718e+04
41.675622.18018.3515e+04
51.599234.93981.3235e+05
Analysis for Stage 5 is finished! ╔═══╗───────╔╗──────╔╗ ║╔══╝───────║║──────║║ ║╚══╦╦═╗╔╦══╣╚═╦══╦═╝║ ║╔══╬╣╔╗╬╣══╣╔╗║║═╣╔╗║ ║║──║║║║║╠══║║║║║═╣╚╝║ ╚╝──╚╩╝╚╩╩══╩╝╚╩══╩══╝ ======================================================================================================================================================
PCA reconstruction MAPE is 0.35238%
PCA-PCE test MAPE is 3.3677%
PCA reconstruction MAPE is 0%
PCA-PCE test MAPE is 0.72036%
Warning: cannot fit gumbel distribution to the given data
Warning: cannot fit gumbelmin distribution to the given data
Starting AIES... |##############################| 100.00% Finished AIES!
Warning: When plotting RSM-response surface for outputfield 1, no observation for this stage!!
Print_MAP = 6×3 table
 K_{0}OCRG_{0}
11.350030.99851.0594e+05
21.373543.65851.0414e+05
31.404120.43716.7718e+04
41.675622.18018.3515e+04
51.599234.93981.3235e+05
61.564132.92168.4291e+04
Analysis for Stage 6 is finished! ╔═══╗───────╔╗──────╔╗ ║╔══╝───────║║──────║║ ║╚══╦╦═╗╔╦══╣╚═╦══╦═╝║ ║╔══╬╣╔╗╬╣══╣╔╗║║═╣╔╗║ ║║──║║║║║╠══║║║║║═╣╚╝║ ╚╝──╚╩╝╚╩╩══╩╝╚╩══╩══╝ ======================================================================================================================================================
PCA reconstruction MAPE is 1.2653%
PCA-PCE test MAPE is 1.4108%
PCA reconstruction MAPE is 0%
PCA-PCE test MAPE is 0.26871%
Warning: Unable to compute a covariance matrix because the computed Hessian matrix is not positive definite.
Warning: cannot fit student distribution to the given data
Warning: cannot fit gumbel distribution to the given data
Warning: cannot fit gumbelmin distribution to the given data
Starting AIES... |##############################| 100.00% Finished AIES!
Warning: When plotting RSM-response surface for outputfield 1, no observation for this stage!!
Print_MAP = 7×3 table
 K_{0}OCRG_{0}
11.350030.99851.0594e+05
21.373543.65851.0414e+05
31.404120.43716.7718e+04
41.675622.18018.3515e+04
51.599234.93981.3235e+05
61.564132.92168.4291e+04
71.528539.84637.5785e+04
Analysis for Stage 7 is finished! ╔═══╗───────╔╗──────╔╗ ║╔══╝───────║║──────║║ ║╚══╦╦═╗╔╦══╣╚═╦══╦═╝║ ║╔══╬╣╔╗╬╣══╣╔╗║║═╣╔╗║ ║║──║║║║║╠══║║║║║═╣╚╝║ ╚╝──╚╩╝╚╩╩══╩╝╚╩══╩══╝ ======================================================================================================================================================
PCA reconstruction MAPE is 0.18715%
PCA-PCE test MAPE is 0.52149%
PCA reconstruction MAPE is 0%
PCA-PCE test MAPE is 0.35515%
Warning: Unable to compute a covariance matrix because the computed Hessian matrix is not positive definite.
Warning: cannot fit gumbel distribution to the given data
Warning: cannot fit gumbelmin distribution to the given data
Starting AIES... |##############################| 100.00% Finished AIES!
Warning: When plotting RSM-response surface for outputfield 1, no observation for this stage!!
Print_MAP = 8×3 table
 K_{0}OCRG_{0}
11.350030.99851.0594e+05
21.373543.65851.0414e+05
31.404120.43716.7718e+04
41.675622.18018.3515e+04
51.599234.93981.3235e+05
61.564132.92168.4291e+04
71.528539.84637.5785e+04
81.467831.54737.7362e+04
Analysis for Stage 8 is finished! ╔═══╗───────╔╗──────╔╗ ║╔══╝───────║║──────║║ ║╚══╦╦═╗╔╦══╣╚═╦══╦═╝║ ║╔══╬╣╔╗╬╣══╣╔╗║║═╣╔╗║ ║║──║║║║║╠══║║║║║═╣╚╝║ ╚╝──╚╩╝╚╩╩══╩╝╚╩══╩══╝ ======================================================================================================================================================
PCA reconstruction MAPE is 0.29337%
PCA-PCE test MAPE is 0.72426%
PCA reconstruction MAPE is 3.6106e-16%
PCA-PCE test MAPE is 0.24529%
Warning: cannot fit gumbel distribution to the given data
Warning: cannot fit gumbelmin distribution to the given data
Starting AIES... |##############################| 100.00% Finished AIES!
Print_MAP = 9×3 table
 K_{0}OCRG_{0}
11.350030.99851.0594e+05
21.373543.65851.0414e+05
31.404120.43716.7718e+04
41.675622.18018.3515e+04
51.599234.93981.3235e+05
61.564132.92168.4291e+04
71.528539.84637.5785e+04
81.467831.54737.7362e+04
91.480637.89861.2207e+05
Analysis for Stage 9 is finished! ╔═══╗───────╔╗──────╔╗ ║╔══╝───────║║──────║║ ║╚══╦╦═╗╔╦══╣╚═╦══╦═╝║ ║╔══╬╣╔╗╬╣══╣╔╗║║═╣╔╗║ ║║──║║║║║╠══║║║║║═╣╚╝║ ╚╝──╚╩╝╚╩╩══╩╝╚╩══╩══╝ ======================================================================================================================================================